Exact BCS stochastic schemes for a time dependent many-body fermionic system 
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The exact quantum state evolution of a fermionic gas with binary interactions is obtained as 
' the stochastic average of BCS-state trajectories. We find the most general Ito stochastic equations 

£\j t which reproduce exactly the dynamics of the system and we obtain some conditions to minimize the 

stochastic spreading of the trajectories in the Hilbert space. The relation between the optimized 
equations and mean-field equations is analyzed. The method is applied to a simple two-site model. 
\ The simulations display effects that cannot be obtained in the mean-field approximation. 

\Q '■ 

I. INTRODUCTION 

The numerical solution of a time dependent many-body problem is a formidable task for large quantum systems. 
For example, for a system of an arbitrary number of fermions with M possible modes of the quantum field, e.g. in a 
lattice model, the dimension of the Hilbert space is 2 M so that both the computer time and the memory requirements 
scale exponentially with the number of modes and becomes rapidly intractable when M increases. A similar situation 
occurs also in classical mechanics, when we describe the dynamics as the evolution of a probability distribution. 
To circumvent the problem on the memory requirement in classical physics, an approach is not to solve numerically 
| the equation of motion of the probability distribution, but to solve the statistical evolution of the variables and 
to evaluate the mean value of some quantities over a finite number of realizations. Such a Monte Carlo approach 
can be used also in quantum mechanics, the most famous example being the Path Integral Monte Carlo based on 
the Feynman's path integral formulation However, apart from notable exceptions (such as the imaginary time 
evolution of bosons with real Hamiltonians, or the imaginary time evolution of some models of fermions [2|), for a 
general quantum problem, the known Quantum Monte Carlo methods do not solve the computer time problem, which 
^ ■ remains exponentially long because an exponentially large number of Monte Carlo realizations is usually required 
For fermions, this problem is the celebrated sign problem, which has been the subject of many efforts, both for 
real-time and imaginary-time simulations pi H 11 M. 1st IsL ITol fill 111 IT^ . 
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In Monte Carlo techniques with path integral, the randomly generated states are generally mutually orthogonal. 
In the original formulation of Feynman, they are the eigenvectors of the particle coordinates. As an alternative, we 
can use to randomly explore an over-complete set of states. Since the dynamics of degenerate bosonic gas with weak 
interactions is approximatively described by the evolution of a Hartree-Fock state, it can be convenient to evaluate 
the exact dynamics using a superposition of paths of Hartree-Fock states . A similar approach has been used for 
fermion systems ^BJ . In Reference 01 the case of bosonic coherent states is also studied. 

The number of random paths necessary to describe the dynamics can be reduced by increasing the number of 
elements of the over-complete set. In the extreme limit where every state of the Hilbert state is an element of the 
explored set, the dynamics can be described with a single, deterministic path: this corresponds to solving directly 
the Schrodinger equation, but this faces again the memory problem. As an intermediate possibility, we can choose a 
£j , set of elements whose single path is a better approximation to the exact solution than the Hartree-Fock ansatz, but 
I* ■ which is still numerically tractable. Attractive interactions in a fermionic gas can lead to the condensation of Cooper 
. £h ! pairs in the superfluid state, as currently investigated experimentally in atomic gases close to a Feshbach resonance 
16]. It is expected that such a superfluid state is reasonably well described by a BCS-state, much better indeed than 
by a Hartree-Fock state. For this reason we study here an exact stochastic approach with BCS states. 

In this article we consider the dynamics in real-time of a system of fermions with binary interactions on a spatial 
lattice. The Hamiltonian is 

fi = ^2,hkic\ci + ]^2v M c\c\cic kl (1) 

kl kl 

h and V are hermitian and real symmetric matrices, respectively, and c^, c\, are Fermi annihilation and creation 
operators. The mode index k labels the spin state <7fc and the lattice node in position rfe. In what follows, we shall 
denote as m s the total number of lattice nodes. 

We wish to obtain the dynamical evolution of the quantum state as the average of stochastic trajectories of BCS 
states. The exact evolution is achieved by averaging an infinite number of stochastic trajectories. The BCS state 
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ansatz that we use is 

|0,7> = ftS(7)|0) ee (lejSw^lo), (2) 

7 being an antisymmetric matrix, involving for spin 1/2 fermions a number of variables 2m s (2m s — l)/2, and fl being 
a multiplicative complex variable. Note that the state in Eq. J2J is in general not normalized. We shall consider the 
case where both 7 and ft are stochastic variables solving Ito stochastic equations 

In Section [H] we find the necessary and sufficient conditions on the stochastic equations in order to have an exact 
description of the dynamics. These constraints do not fix univocally the stochastic scheme, thus we shall use this 
freedom to reduce the statistical spreading of the trajectories. In Section ITTT1 we construct explicit stochastic schemes. 
The growth rate of the spreading is evaluated and an upper limit for the statistical error on the observables is 
established, which shows that the statistical uncertainty is finite at every finite time. In Section llVl the stochastic 
approach is illustrated on a two-site model. 



II. STOCHASTIC EQUATIONS 



A. Conditions for the stochastic evolution to be exact 



We want to evaluate exactly the quantum state evolution using a superposition of the BCS states |ft, 7) = 0| 7 ) 
with a stochastic evolution of 7 and ft. For an infinitesimal variation of 7 and ft we calculate the variation of the 
ansatz by expanding |ft + Aft, 7 + A 7 ) in powers of Aft and A7: we have from Eq. I|A4|I that 



A|n )7 > 



Ail 1 
7T + 2 



ijkl 



|n.7>, 



(3) 



On the other side, the Hamiltonian evolution during At of the state equal to |ft,7) at time t is given to first order in 
At by Schrodinger's equation: 



iHAt|Q,7> 



\ E V ^kljic\c]clc\ At - iE ( -^Ulii + E hik ^W ) At 

ijkl ij \ k / 



(4) 



where we used Eq. (|A5|) to express c|ft, 7 ) in terms of 7} and where we took h = 1. If 7 and ft satisfy a 
deterministic equation, it is obvious that the first term of the right hand side of Eq. |0J does not in general coincide 
with the third term of the right hand side of Eq. (©. As we shall prove, Eq.Q and Eq.© can become equal when we 
consider stochastic equations and we average Eq. over every possible realization of the stochastic variation during 
At: 



A|ft, 7 ) = -i«Ai|n,7). 

Since S^) is invertible, we have to find a stochastic equation for ft and 7 that satisfies the equality 



(5) 



ij ijkl ij 

\ E v ^ikijic\cj d l4 At - ! E ( 2 Vijlij + hik7kj ) 4 £ l At 

ijkl ij \ k / 



|o> 



|0). 



Equation (jfjj) is equivalent to the following ones, 



Aft 

7T 



= 



A7y = -^AQA'jij - iVijjijAt - iE h lk j kj At + i^ h jk j ki At 



E 

permutation of ijkl 



-AjtjAju - -VijjikjjiAt 



(6) 

(7) 
(8) 

(9) 
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where (— l) p is the signature of the permutation. The first equation implies that the deterministic term of n is zero. 
The second equation gives the deterministic term for 7, the last one gives a condition for the noise term of 7. This 
last condition can be written explicitly 



A7yA 7fci + A ljk A la + A 7ife A7y + iAi(7« + V kj + V u + VuHklij 
+iAt(V ik + V kj + Va + V 3 i) ll]lk i + iAt{V lk + Vik + Vy + Vy)7jfc7ii = (10) 

Note that this equation is automatically fulfilled when two of the four indices ijkl are equal. 



B. Growth of the statistical error 



To estimate the statistical error of the method, we consider the growth rat e of the mean sq uared distance between 
the true state of the system and a single realization of the stochastic ansatz, A||^>) — |S1,7)| 2 = AM, where 

M=(0, 7 |0,7). (11) 

This will allow to prove that the statistical error remains finite at all finite evolution times and this will provide a 
strategy to identify optimal stochastic schemes in trying to minimize the growth rate AM /M . 
To the first order in At, 



AM = (A<n )7 |)|fi,7> + <n, 7 |(A|fi,7» + (A(0, 7 |)(A|n > 7)). 



(12) 



In the right hand side, the sum of the first two terms gives exactly zero, since by construction A|S1, 7) = — iAtH\£l, j)/h. 
In the last term, we can replace A|f2,7) by its stochastic component: 



A|Q,7> 



stoch 



A|0, 7 )-A|n,7) 
9. 



|n.7>. 



(13) 



where we used the fact that Aft is purely stochastic and where A7j Stoch is the stochastic part of A7y . Equation 1|12|) 



leads to 



AM ||A|f7,7) stoch || 2 



M 



(n, 7 |n )7 ) 



(14) 



which can be evaluated using Wick's theorem: 



AM 
~M 



Ail 



|E A ^ A 7fc*(e4>(^4>- 



(15) 



ijkl 



where the expectation value is taken in the ansatz, (...) = (^,7| • ■ ■ \Q,"f)/M. A clear step for the minimization of 
the error growth is to choose Aft in order to set to zero the first term in the right hand side of the above expression: 



Afl 

~~n~ 



kl 



We shall always choose An in this way in what follows. It is then easy to show [see Eq. I)13|l] that 



AM = AM, 



(16) 



(17) 



i.e., the stochastic terms of AM are exactly zero. Furthermore the deterministic part of A7 is now slaved to the 
stochastic part of A7, according to Eq.(|5J): the only increments that remain to be specified to fully determine the 
stochastic scheme are A7|j och , and this we shall do in the next section. 
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III. EXPLICIT EXACT STOCHASTIC SCHEMES 
A. Our solution for an arbitrary interaction potential 

This most general solution relies on the following ansatz for the stochastic increment: 

A 7 g° ch = (Afi + Afj)jij (18) 

where the noise terms A/& are independent of 7^ . Inserting this ansatz in the validity condition Eq. (|10|l , we find that 
if the noise terms have the following correlation function, 



A/, Af 3 = -iVyAi, (19) 

this validity condition is satisfied |2lj| . Since the matrix V%j is real symmetric it can be diagonalized; a noise having 
this correlation function may then be explicitly constructed using the corresponding eigenbasis. 
In the specific case of a discrete 8 interaction potential between two opposite spin components: 

Vij = V^ rt<rj S au - a} (20) 

where r, and Ui are the lattice position and the spin component of the mode of index i, the following explicit noise 
may be used: 

fi = (-iV ) 1/2 [A£ ri ^ j>T + A£M,l] (21) 

where the A£ r 's are statistically independent complex Gaussian noises of variance At. For this specific noise imple- 
mentation, f*fj = \Vo\AtSij, so that the growth rate of the statistical error can be expressed by the simple formula 

M 

^ = \V \AtJ2(Zki)(&k) < hv \Atm 8 (22) 
fc 

where m s is the number of lattice nodes and where we used the Eqs. 1JA8IA11IIA9I) . AM/ (MAt) has a constant as an 
upper bound, thus the norm squared of |fi, 7) is bounded at every time by exp(|Vb|m s i/2) times its initial value. A 
similar bound was derived in the stochastic Hartree-Fock scheme in |15| . with a larger exponent Il8l. Consequently, 
the Monte Carlo statistical variance of an observable O is finite when Tr[(9 2 ] is finite [see Ref. [l||, Eq. (30)]. 
To be complete we also give the corresponding deterministic part of the evolution of 7: 

A-yij/At = -iVijfij - i^2{h lk j kj + hjklik) - i^2{V ik + v ik){c\c k )^ij. (23) 

k k 

We note that the last sum over k in the righthand side is simply the Hartree mean-field term. 

B. Case of an off-diagonal 7 and a delta interaction 

We now restrict to useful limiting cases where the one-body Hamiltonian h is spin-diagonal, the interaction potential 
is an on-site discrete <5 between two distinct spin components, as defined in Ea. (|20|l . and where the matrix 7 in the 
ansatz has initially zero matrix elements between identical spin components. We have then identified exact stochastic 
schemes that preserve at any time this block off-diagonal structure of 7. 

1. The solution that we have found with minimal error growth 

A general strategy to find the 'best' stochastic scheme among a very large number of possibilities is to try to 
minimize the growth rate of the statistical error Eq. I|14|) . Whereas this program is easily fulfilled for bosons [l4( it 
seems to be more difficult to achieve for the stochastic BCS ansatz. Here we report the solution that we have found 
with minimal error growth. It is possible that a better solution exists. The stochastic increment is given by 

A 7 f- i h rj = (iy ) 1 /2 (A ^ i - A£ rj ) 7Tn , ir . (24) 
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where the A£ r are m s independent real Gaussian noises of variance At. The corresponding growth rate of the 
statistical error is exactly given by 



AM 



At|F | ^2 [^\Ar){c^ r c\ r ) + (c\ r c lr )(c lr c\ r ) - 2(c| r cj r )(c ir c Tr ) 



(25) 



which is indeed smaller than the general result Eq. l|22|l because of the occurrence of a negative term involving 
anomalous averages. In this scheme the deterministic part of the evolution of 7 is given by 

(26) 



2. The solution that we have found with minimal memory requirement 



In the case when the one-body Hamiltonian h is totally spin independent and when the off-diagonal block 7jr;,ir- is 
initially a symmetric or antisymmetric matrix (under the exchange of rj and rj), we have found a stochastic scheme 
which preserves this symmetry property at all times, which allows to save a factor of two on the memory requirement; 



= *(^o) 1/2 (A6, + A^)7Tra^ 



(27) 



where the A£ r are m s independent real Gaussian noises of variance At. In this case, the growth rate of the statistical 
error is 



AM 



A *N [(c\Ar)(c lr c] r ) + (c\ r C lr )(c lr c\ r ) +2(4 r 4r)(C|rC Tr )] < At|V&|» 



(28) 



which is larger than for the two previous schemes. In this scheme the deterministic part of the evolution of 7 is given 
by 

(29) 



C. Link with the mean-field approximation 



In the three explicit stochastic schemes given in this section, the deterministic part of the evolution for 7^ did not 
coincide with the mean-field evolution. This in contrast with the optimized stochastic Hartree-Fock schemes obtained 
for bosons [l4| and for fermions 0]. For a general stochastic BCS ansatz with the optimizing choice Ea. l|16|l . we found 
the following relation between the deterministic evolution of 7^ and its mean field evolution (given in Appendix : 

A7- - A 7j I f an fi ° ld = E A 7ifcA 7i *<44> (30) 

kl 

by inserting the expression Eq. (|16|l of Ail into Eq. (JSJ and by using Eq. I|1U|I to eliminate AjijAjki- This shows that 
finding a stochastic scheme where the deterministic and mean field evolutions coincide is not straightforward. 



IV. STOCHASTIC APPROACH FOR A TWO-SITE SYSTEM 



In order to illustrate the method, we apply it to a simple system with two sites, corresponding to the Hamiltonian 



H = 



5 £[4. 



C<x2 + cl 2 C<rl 



V 



(31) 



where the spin index a takes the values f and j. There is no interparticle interaction when the two particles are in 
different wells. A physical system that may be described by this model is a set of two Fermi particles in a double-well 
potential. 
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FIG. 1: Mean value of the population in the state 1 | as a function of time. The solid and dashed lines are evaluated using 
the schemes of £1111 B H and illllB 21 respectively. The number of realizations is 10 5 and V = 0.2. The dashed-dotted line is the 
BCS mean-field prediction. The widths of the error bars are the standard deviations. 

At the initial time, we choose a BCS state with the elements of 7 equal to zero, apart from 7if ,14 = — 711, it = 7o = 2. 
The state is a superposition of the vacuum and the state with two atoms in the site 1. The direct numerical solution 
of the dynamics is obtained writing the Hamiltonian in the basis of the Fock states of the operators Ci tS and cj a . 
The integration is simplified by the fact that the interaction cannot flip the spin and the amplitude of some states 
remains zero. In Fig. ^ we report the mean value of the population in the state f 1 as a function of time. We have set 
Vq = 0.2. The dashed-dotted line is evaluated with the mean-field equations (as given in Appendix |B|) . the dashed 
line is the direct numerical solution and the solid line is the stochastic solution. The widths of the error bars are 
the standard deviations. We have used 10 5 realizations with the scheme of 3III B II In the mean-field approximation, 
the evolution has a damped oscillation with a revival for t > 30. The collapse and revival of the oscillations of the 
exact solution occur with a shorter time scale. The stochastic approach is able to display very well this behavior. In 
Fig. [21 we report (£1,7|£1,7) as a function of time. The solid dashed and dotted lines are evaluated using the schemes 
of qill B II OTA1 and 3111 B 21 respectively. The dashed-dotted line is the upper bound of Eq. (|25|l for the first two 
schemes. As expected, the optimized scheme has a smaller spreading. 

Note that the growth rate of M is zero at the initial time for the scheme of 3111 B II (see inset of Fig. 0), because 
of the presence of the last term in Eq. (|25(l . which cancels the other contributions in the initial state considered 
here of particles localized on a site. The spreading of the trajectories grows exponentially and it increases for larger 
interparticle interactions. We have done similar calculations for a stronger interaction, e.g. for V — 2; the stochastic 
method agrees with the direct numerical solution, with a higher growth rate of the statistical error for increasing V, 
as expected. 



In this article we have shown that the state evolution of a fermionic gas with binary interactions can be obtained 
in an exact way as the average of stochastic trajectories of BCS states. We have derived the general Ito stochastic 
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FIG. 2: Average of M — (fi, 7|fi, 7) over the stochastic realisations, as a function of time, evaluated using the scheme of 3111 B II 
(solid line), of 3111 Al (dashed line) and of 3111 B 21 f dotted line). The dashed-dotted line is the upper bound of Eq. 1251 for the 
first two schemes. At the initial time the growth of M is zero for the first scheme (see inset). 



equations which give the exact evolution of the system and we have found a condition on some parameters of these 
equations to reduce the statistical spreading of the trajectories in the Hilbcrt space. The upper bound that we have 
found on the spreading for a particular scheme is similar to the one obtained for the Hartree-Fock ansatz in , with 
a smaller value. We have illustrated the method on a two-site model and we have shown that the quantum effects, 
which cannot be obtained with a mean-field approximation, are displayed by the results of the stochastic approach. 
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APPENDIX A: SOME PROPERTIES OF BCS STATES 



First we prove that the BCS states form a complete family for the states with an even number of atoms. The set 
of states with a definite number of atoms in each mode constitutes an orthonormal basis of the Hilbert space. It is 
sufficient to show that each element of this set is equal to a superposition of BCS states. An element has the following 
form 

N p 

\{k n ,l n }) = (1111)10), (Al) 
n=l 
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where we have grouped the atoms in pairs. The n-th pair has the atoms in the k n and l n modes. N p is the number 
of pairs. k n ^ k m ^ l n =/= l m when n ^ m, whereas k n ^ l m for every n and m. It is easy to prove that 

\{k n ,l n })=K #! / # 2 ... / d^ Np e-^^^ +cx ^< & 0\0) } (A2) 
Jo Jo Jo 

where if is a normalization constant. Thus, the BCS states form a complete family. Actually it is over-complete. 
It is possible to demonstrate that (see Section 2.2 of Ref. ,20j) [2^ | 



and 



From Eq. (|A3|) we have 



Using Eqs. l|A4IA6jl we find that 



M = (n,7|n,7> = iradct[l + 7 t 7 ] 1 / 2 (A3) 



^L\n, 7 ) = c{c}\n,j), (A4) 



d 



: |a7)=^^4l fi .7)- (A5) 



(n,7|c t 4|n, 7 ) 
<n,7|n,7> 

(n, 7 |ctct|n, 7 ) 
(n,7|n,7> 

(^, 7 |CfcC;|»,7) 

(0,710,7) 

fttAtl 



iQ- = -[ 7 t(i+77t)- 1 ] ii M. (A6) 



(fi,7|0,7) 



= [(1 + ^7)-% (A7) 
= [7(l + 7 t 7)" 1 7 t ]ft (A8) 
= -[7(l + 7 t 7)" 1 ]fci (A9) 

(A10) 



We note that Equation (|A8|) can be written in various forms using the matrix identities 

7*7(1 + 7 t 7)~ 1 = 7 t (l + 77 t )~ 1 7 = [7(1 + 7 t 7) _1 7 t ] T (AH) 
where A T is the transpose of matrix A. 



APPENDIX B: MEAN-FIELD EQUATIONS 



Using the results of section 9.9b of |2fl|, we obtain the following equations of motion for 7 in the mean-field 
approximation: 



la 



■ i Vu {ckCiYlkilji + Wjj (cjCj). 

(Bl) 
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